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A fermionic operator circuit is a product of fermionic operators of usually different and partially overlap- 
ping support. Further elements of fermionic operator circuits (FOCs) are partial traces and partial projections. 
The presented framework allows for the introduction of fermionic versions of known qudit operator circuits 
(QUOC), important for the simulation of strongly correlated <i-dimensional systems: The multiscale entangle- 
ment renormalization ansatze (MERA), tree tensor networks (TTN), projected entangled pair states (PEPS), or 
their infinite-size versions (iPEPS etc.). After the definition of a FOC, we present a method to contract it with 
the same computation and memory requirements as a corresponding QUOC, for which all fermionic operators 
are replaced by qudit operators of identical dimension. A given scheme for contracting the QUOC relates to 
an analogous scheme for the corresponding fermionic circuit, where additional marginal computational costs 
arise only from reordering of modes for operators occurring in intermediate stages of the contraction. Our result 
hence generalizes efficient schemes for the simulation of d-dimensional spin systems, as MERA, TTN, or PEPS 
to the fermionic case. 

PACS numbers: 03.67.-a 02.70.-c, 71.10.Fd, 



I. INTRODUCTION 

Strongly correlated quantum lattice models pose some of 
the most intriguing physical questions and technical chal- 
lenges, due to the fact that the number of degrees of free- 
dom increases exponentially with the system size. Classify- 
ing the intricacy of calculating ground state energies of such 
systems has become a vivid branch of complexity theory 03- 
|3| . Especially for the analysis of ground state properties in 
one-dimensional systems, the density -matrix renormalization- 
group (DMRG) [4, 5 | provides a numerical approach that is 
often extraordinarily accurate. It works by variational opti- 
mization of a suitable class of states, so-called matrix product 
states 1 6-8 1 . For two- and three-dimensional systems, quan- 
tum Monte-Carlo methods (e.g., positive-definite path integral 
||9] [TO) or stochastic series expansion ifTTI representation) are 
extremely successful for bosonic and unfrustrated spin mod- 
els, but are bothered by the sign problem ifTOlfTZl for some in- 
teresting frustrated spin and fermionic models, including the 
notorious Fermi-Hubbard model 

H = - Yl (fLfja + h - c -) + U Y1 n^hn - fiJ2n ia 

which is a candidate for the description of the essential 
physics of high-temperature superconductivity. Recently, new 
tools such as the diagrammatic Monte Carlo method have 
been developed fT3l[T4"l . which have a less severe sign prob- 
lem and have, e.g., been demonstrated to give precise results 
for the repulsive Fermi-Hubbard model in the (correlated) 
Fermi liquid regime IT31 . 

In a complementary development, generalizations of 
DMRG ideas to higher dimensions have been put forward. To 
this purpose, first, one needs to give an ansatz for the many- 
particle state for which the number of degrees of freedom does 
only scale polynomial with system size but is (hopefully) still 
appropriate to describe, e.g., the ground states of the higher- 
dimensional system. Second, a way of efficiently evaluating 



interesting local observables or correlators with respect to the 
ansatz states needs to be identified. Third, a corresponding al- 
gorithm to determine or approximate the ground state within 
the ansatz class on a classical computer needs to be worked 
out. Focusing first on spin (or equivalently qudit) lattices, sev- 
eral suggestions have been put forward, such as tensor prod- 
uct ansatze or projected entangled pair states (PEPS) II 1 6442 1 1 , 
tree tensor networks (TTN) 1221 . or multiscale entanglement 
renormalization ansatze (MERA) 112344271 . 

In this article we address the question of how higher- 
dimensional fermionic systems can be studied via ansatz 
states. If one maps the system to a spin model by express- 
ing states and operators in the occupation number representa- 
tion with respect to a fixed ordering of the modes, inevitably 
long-range (0(L d_1 ), where L is the linear size of the d- 
dimensional lattice) interaction terms occur, rendering sim- 
ulation unfeasible: The spin representation of a term fjfk, 
j < k, under the Jordan-Wigner transformation |28| is for 
instance 

o-J® (g) <rf®a+, 

j<l<k 

containing a so-called Jordan-Wigner string. 

Accompanied by first numeric results, very recently, 
fermionic generalizations of MERA states were suggested in 
Refs. lHUiol and for PEPS in Ref. EJJ. Specifically, in Ref. 
OUll . also an algorithm for fermionic MERA is given by dy- 
namical reordering. It exploits the possibility to change the 
ordering of the fermionic modes during the algorithm to con- 
fine all occurring Jordan-Wigner strings to a sublattice of fi- 
nite extent, the causal cone of, e.g., a local observable in the 
MERA. Going beyond that result, here, we pose the question 
whether a given general circuit of fermionic operators (FOC, 
examples in Fig.[JJ) can be contracted with the same efficiency 
as a corresponding circuit of qudit operators (QUOC). This is 
answered in the affirmative for the case where each operator 
in the FOC is parity-symmetric (either fermion number par- 
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Figure 1 : (a) The graphical representation of a FOC as a directed graph. The nodes represent fermionic operators. The arcs (directed edges) 
represent (partial) multiplications, partial traces, and open indices. Each arc is labeled by the set of modes it corresponds to. The operator 
corresponding to a certain vertex maps from the modes of all incoming arcs to the modes of the outgoing arcs. In the example, the arc "e" 
corresponds to a partial multiplication, arc "p" to a partial trace, and arcs "a" and "6" to open incoming and outgoing indices, respectively. The 
node at the top corresponds to a ket vector from J- m un and the node at the bottom left to a bra vector (element of the dual of -F/uiuj)- As a 
whole, the circuit is a fermionic operator mapping from T a to Th- (b) A FOC for the calculation of the expectation value of a local observable 
(square in the center) with respect to a MERA state with two renormalization steps. The hatched flat rectangles represent isometries which 
correspond to a coarse graining step in a (realspace) renormalization procedure. The other rectangles represent unitaries that are supposed to 
reduce entanglement of adjacent blocks before a coarse graining step. The circuit contains only those unitaries and isometries of the MERA 
that lie inside the so-called causal cone of the observable; all others cancel out. (c) FOC for a tree tensor network (TTN) state, here for a 
genuine tree system, the Bethe lattice with coordination number z — 3. To have the value of a FOC well-defined, one needs to specify an 
ordering among the operators, assigning to each operator a number r = 1, 2, .... In example (a), we arbitrarily chose r to increase from the 
bottom to the top. In example (b), a natural ordering, motivated by the picture of subsequent renormalization steps, is also directed from the 
bottom to the top; as we will explain later, the ordering inside one layer is irrelevant, as the contained isometries are all parity-preserving and 
operate on disjoint sets of modes. Analogously in example (c), we can choose r to increase in radial direction, starting from the central node. 



ity preserving or changing): We show constructively that the 
elementary contraction operations for such a FOC can be ex- 
ecuted in an arbitrary sequence and give a detailed account of 
the algorithm. As compared to the requirements for the con- 
traction of a certain QUOC with a given contraction scheme, 
the number of operations and memory requirements for the 
same contraction scheme, applied to a corresponding FOC, 
increase only by a marginal amount. 

This allows to translate the algorithms already developed 
for spin systems (for PEPS, e.g., in Refs. lfT71[T9ll32H34l. for 
MERA, e.g., in Refs. J24] |26] W\ E2)) to the fermionic case 
without loss of computational efficiency. Giving further de- 
tails for the case of PEPS, we argue that application of the 
FOC scheme to fermionic PEPS appears to provide a more 
efficient algorithm than that presented in Ref. ||3T| where a 
mapping to a spin system was employed by choice of a fixed 
mode ordering. 

In Sec. |ll]the idea of the FOC is introduced and it is given 
a proper definition. Rules for the execution of the elementary 
contraction operations for two or one operators are derived 



in Sec. Ill after which the importance of a predefined order 
among the operators constituting the FOC is pointed out in 



Sec. IV It is also explained how this operator order can be 
modified with marginal computational cost, allowing to ef- 
ficiently execute the elementary contractions in an arbitrary 
sequence. The implications on computational efficiency and 
locality considerations are summarized in Sec.[V] Sec. VI 



in- 



troduces further useful operations on FOCs that are employed 
in an efficient contraction algorithm for fermionic (i)PEPS in 
Sec.lVIII The article closes with a short discussion. 



II. FERMIONIC OPERATOR CIRCUIT 
A. General structure 

A fermionic operator circuit (FOC) is a product of (not nec- 
essarily physical, i.e., in general not particle number parity 
preserving) fermionic operators A t : T m T n of in general 
different support, specified by sets of mode labels m, n C C. 
Further elements of FOCs are partial traces and partial projec- 
tions. Each mode label x E C occurs at most twice, once for 
an incoming mode of some operator and, the second time, for 
an outgoing mode of the same or another operator. This means 
for graphical representations of FOCs as graphs, where each 
vertex corresponds to one operator Ai, that each arc (directed 
edge) of the graph carries a set of unique mode labels. As 
explained in Sec. |II B| this convention allows for a convenient 
definition of the FOC such that it has a well-defined value. 

Prominent examples of FOCs are fermionic versions of 
known qudit operator circuits (QUOC), important for the sim- 
ulation of strongly correlated rf-dimensional systems: mul- 
tiscale entanglement renormalization ansatze (MERA) [23 II 
and tree tensor networks (TTN) |22|; Fig. [T] As we show in 
Sec. |VII| also the fermionic variants of tensor product ansatze 
or projected entangled pair states (PEPS) lfl6HT9l are covered 
in the FOC framework; Fig. [8] For a MERA, a possible choice 
for mode labels are the renormalization step t combined with 
a site label from the corresponding lattice. 

For numerical purposes, each fermionic operator A : 
J- m — > J- n of the circuit is stored in an occupation num- 
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ber representation with respect to certain orderings m and n 
of the sets of modes m,n C C. We consider such order- 
ings as bijective enumerations m : {1, . . . , \m\} — > m and 
n : {1, . . . , \n\} — > n of the sets, where \m\ denotes the num- 
ber of elements in m. We may also treat such enumerations as 
vectors. For a chosen ordering n of the modes in n, we denote 
the basis states of the Fock space T n by 

\n) n =K,...,n w ) n :=(4) ni ...(^ [B| r'-'|0)n, tt) 

where |0)„ labels the vacuum state of the Fock space T n and 
fi are the corresponding anticommuting ladder operators with 
{fii fj} — The operator A can hence be stored as the 
complex 2 1 "I x 2l m l matrix 



J n , m (A) = ^2 \n) n (n\A\m) m (m|. 



(2) 




This is an occupation number representation or Jordan- 
Wigner transform ll28l of the operator A. Of course it is 
also possible to restrict (for each set of modes) to a reduced 
basis. The only information about the basis states actually 
needed is their particle number parity; see Sec. |III A The 
states occurring in |2} are elements of different Hilbert spaces: 



e T n , \m) e Si m i, and \n) e B\ n \, 



where Bi n i denotes the |n|-qubit Hilbert space 



B 



\n\ 



2\<g>|n| 



(3) 



A similar approach can be used for anyonic systems 



Figure 2: (Color online) The operator order goes from the bot- 
tom to the top. Left: Example for an operator circuit on a lattice 
L — mUnUpUq. It corresponds to the expression Tr mU q( n(0|-44 ■ 
... • Ai\0) p ), cmp. Eq. |4f, where, e.g., Ai = di ® ld p u q - 
For convenience, we require in the definition of fermionic opera- 
tor circuits, Sec. |IIB| that each mode occurs at most twice, once 
as an incoming mode and once as an outgoing mode of some op- 
erators. This can be achieved by a relabeling of the modes, yield- 
ing the FOC Ai o . . . o Ai (right). This does not change the ma- 
trix elements of the operators and the FOC. One has for example 



m< 2 )en< 2 > 



(m'n'\Ai\mn) 



„{m' n' \ai\mn) 



men > 



where m, n, m , n {1> are orderings of the sets of modes m, n, rrii, 
and Hi. Here, with the relabeling, also the partial projections for 
operators a,2 and az have been executed. 



B. Definition of a FOC 



C. Remarks on the definition 



A fermionic operator circuit is specified by a set of 
fermionic operators {Ai : T mi — > J- ni }, where each mode 
label occurs at most twice, once as an incoming mode of an 
operator Aj and once as an outgoing mode of an operator Aj . 
Mode labels which occur two times in this fashion imply a 
(partial) multiplication, Fig. |4^, or (partial) trace, Fig. [4j}, of 
the corresponding operators with respect to that set of modes. 
Both operations together define a general contraction of two 
operators, namely contraction of some outgoing modes of A 
with some incoming modes of B and, simultaneously, of some 
incoming modes of A with some outgoing modes of B; see 
Fig. |4j;. Mode labels, which occur only once, correspond to 
modes that the FOC as a whole maps from or maps to. 

To have the value of a FOC well-defined, one needs to spec- 
ify an ordering of the contained operators {Ai}. The value of 
the FOC is then defined by the one resulting from doing the 
contractions in the order A^ o . . . o A 2 ° Ax, where "B o A" 
denotes the contraction of all common modes of the operators 
A and B; see Fig.|4j;. 



As discussed in Sec. 



IV 



this operation 

is associative but in general not commutative, Bo A ^ A o B. 



In Sec. [Ill] as for the partial contraction operation, we will 
also give a rule for a partial projection of some modes to basis 
states (i.e., {hi} eigenstates). This is actually already covered 
by the contraction operation but perhaps useful to have explic- 
itly, as such projections are frequently used in considerations 
on operator circuits. 

Note that the operators {Ai} are not assumed to be from the 
so-called algebra of physical operators - i.e., particle number 
parity preserving. This is for example useful when calculating 
correlators of the form (fj fj) with respect to MERA or TTN 
states. In such a calculation, the operators fj and fj become 
(clearly not parity preserving) elements of a FOC. 



However, it will be explained in Sec. IV that in order to be 



able to do the contraction of the FOC in an arbitrary sequence 
(necessary to get optimum numerical efficiency), i.e., to be 
able to deviate from the order An ° • • • ° A% o Ax, it is in 
general necessary that each Aj is either parity preserving or 
parity changing. 

That mode labels are required to be unique is not a limita- 
tion. Consider for example an operator circuit that is defined 
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on a lattice L and does not have that property, 



Tr t ( o (0|^AT 



A 2 -A 1 \0) l 



(4) 



Here t C L denotes a subset of modes that are traced out, and 
i and o C L denote subsets of modes that are projected out; 
t n (i U o) = 0. The circuit hence maps from J^LKttui) to 
•FjA(tUoV Each operator .Aj acts nontrivially on a subset of 
the modes: Ai = ai® Id^x^ with fij : Ti i — > where li C 
L. Now, relabeling of the modes to make the modes unique 
as depicted in Fig. [2] does of course not change the matrix 
elements of the FOC. It yields a proper FOC A^o. . .0A20A1, 
where each operator Ai has the same matrix elements as the 
corresponding hi (partial projections onto the vacuum can be 
executed in the same step, as in our example, or introduced 
as separate elements of the FOC). The contraction rules in 
Sec. [Ill] are constructed such that this FOC and Q have the 
same matrix elements, i.e., are related by a trivial relabeling 
of incoming and outgoing modes. 



The number of particles in a basis state \m) m is denominated 
by 



m := > m 



(6) 



The parity of the basis state is (— l) m . 

Whenever we refer to Fock spaces for unions of sets of 
modes, as in F m \j n , it is implied that those sets of modes are 
disjoint, i.e., m n n = in that case. 

With B •„ A, a partial multiplication is denoted. Only the 
outgoing modes n of A are contracted with the correspond- 
ing same incoming modes n of B. Correspondingly Tr r B 
denotes a partial trace, the contraction of incoming modes r 
with outgoing modes r. By B o A, we denote a (partial) con- 
traction of all common incoming/outgoing modes of A with 
corresponding outgoing/incoming modes of B. 



III. CONTRACTIONS 



D. Rationale behind calculations and derivations 

• The fermionic operators are maps from one Fock space 
of "incoming modes" to another (in general unrelated) 
Fock space of "outgoing modes". In general, they are 
of different dimension. 

• Each arc (directed edge) in a graphical representation of 
a FOC corresponds to a set of unique fermionic modes. 

• Vacuum states are mode specific. Ladder operators of 
other unrelated modes commute with the vacuum state 
for other modes. Take for example n = {1, 2} and n = 
(1,2), then 

fl\n in2 ) n = fMr(fl) n2 \0)n 

= (-l) ni+ ^(/i") Bl (#H0>n -fl 



(-ir + ">m 2 ) n -/ 3 T 



(5) 



The rationale behind this is that if we have an expres- 
sion m un(0\A m A n \0)mun for disjoint sets of modes 
771 and n, and where A m and A n are polynomials in the 
ladder operators of the modes in m and n, respectively, 
we have 

i(0\ m(0\A m A n \0) m |0) 

n rn m n 

(0\A n \0) 

1/ 



E. Notation 



In the following, rules are given for all elementary con- 
traction operations needed during the evaluation of a FOC. 
No non-local Jordan- Wigner transformations occur. The only 
reordering of modes necessary is for incoming or outgoing 
modes of single operators, directly before a partial multiplica- 
tion, trace etc. that they are affected by. 



A. Reordering of modes 

Assume we are given a fermionic operator A : J- m — > T n 
in the occupation number representation J n , m (A). The con- 
traction rules to follow, will pose some preconditions on the 
orderings of modes (to get simple formulae). We need hence 
to be able to derive from J n ,m{A) representations J n >,m'{A) 
with different mode orders. 

All reorderings can be written as sequences of two 
mode swaps. Let us assume that m' = m and that or- 
derings n and n' differ only in modes tv, and r\k (for 



1 < 



J < 



< 



MX 



i.e. 



= n/j and n',, = 



xij. Where in the old representation, \n) corresponds to 
the state |n) n = (/ ni )"^ . . . (/ n .)^ . . . (/ tt J»» . . . | )„ , 
it corresponds in the new representation to |n) n / = 

(/n 1 ) ni ...(/nj^...(/n j )" fc ...|0)„. 

To derive the corresponding transformation on the represen- 
tations of A, note that an operator Sjk that swaps the modes, 
i.e., S jk fjS] k = f k and S jk f k St = f 3 is given by W\ 



We use the Einstein summation convention, i.e., basis state 
labels that occur twice in an expression presuppose summa- 
tion over that basis. 

Basis states for a certain set m C C of \m\ fermionic modes 
and an ordering m of those modes will be denoted by \m) m = 
(F™)t|0) m , where m G {0, l}l m l and 



F r , 



(frc 



(fm 2 ) 2 (/mi) 



s j k = i-f}f j -flh + flf k + fV j - ( 7 ) 

With Sjk\0) n = \0)n > we have hence 

The occupation number representation (Jordan- Wigner trans- 
form) of a term fjf k is a~ ® (®j=j+i of) (8 where the 
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■ — 



tli n 2 n 3 




orderings Ci = I © p, c 2 = m © q is 

C = B- n A 

= \k)t (k\B\n'q) b (n'q\ ■ \np) a (np\A\m) m (m\ 
= (_ 1 )M+(P+9)(«+rV)| fc ^ { (k\B\n'q) 

x q (0\ n (0\(.FgyF?'(F?yF*\0) n \0) p 

x(np\A\m) m (m\ 
= |fep) Cl (fe|B|nq)(np|A|m) C2 (m 9 | 

=: IMdCM^lTngJcaM, (10) 

where C is the representation C = J Cl>C2 (C). In short, the 
transformation rule for the occupation number representations 
reads 



Figure 3: (Color online) To implement contraction schemes for FOCs 
on a computer, we represent every operator A : T m — > T n in an oc- 
cupation number representation J n .m(A). The primitive contraction 
rules, given in Sec. [Ill] pose some preconditions on the orderings of 
modes (to get simple formulae). Hence, before applying those rules, 
it is in general necessary to change, e.g., from J n ,m(A) to a represen- 
tation J„', m i (A) with different mode ordering. In the depicted exam- 
ple, the order of the outgoing modes chan ges fro m n = (ni, ri2, rt3) 
to n' = (ri2, ri3, ni). As explained in Sec. Ill A this requires appli- 



cation of the swap matrix S [Eq. {9}] - in this example two times 



(kp\C\mq) = {-l) m {k\B\nq){np\A\m). 



(11) 



In appendix |Aj an alternative derivation of this rule is given, 
where the support of operators A and B is extended prior to 
the contraction such that there is no need for applying the 
commutation prescription Q. The result is the same. 



C. Partial trace of an operator 



a a denote the Pauli matrices. The swap operator for two con- 
secutive modes is in the relevant subspace 

S '■= J(i,i+l),(i,i+l){Si,i+l) 

= |0,0)(0,0|-|1,1)(1,1| 

+ |0,1)(1,0| + |1,0)(0,1|. (9) 

In practice one may choose to execute all mode reorderings 
by application of corresponding sequences of swap operators 
for consecutive modes; see Fig. [3] 

Swapping of whole sets of modes, e.g., useful when retain- 
ing reduced bases, can be done as well. Consider an operator 
B : F m -> F uUvUxUz given in the representation J n , m (B) 
with n = uffit>ffir.ffi3 where u, 0, p, 3 are orderings for the 
modes in u, v, x, and z. Swapping c and y, is achieved by 
(uxvz\J n i , m (B)\m) = {— \) xv (uvxz \ J n , m (B)\m), where 
n' = u©y © ©3. 



B. Contraction of some outgoing modes of A with the 
corresponding incoming modes of B 

The partial multiplication of two operators is depicted in 
Fig.|4^. Let A : T m — > F n up and B : J"„u 9 — > Fk, i-e., the 
operators' outgoing/incoming supports overlap in the modes 
n. Let m, n, p, and q be orderings for the modes in m, n, p, 
and q. Assuming we have the two operators in representations 
A = J a , m (A) irndB = J {jb (-B) with a = nffip and b = nffiq, 
the resulting operator C := B •„ A : J- m uq J~kup with 



The partial trace of an operator is depicted in Fig. [4j}. Let 
A : J mUr — > J~ nUr , i.e., the operator's outgoing and incom- 
ing supports overlap in the modes r. Such operators can al- 
ways be decomposed in the form 

i = + (12) 

where A + is the particle number parity preserving and A- the 
parity changing component, i.e., 

(_l)*»+*rXt = ±A ± (-lf™+^ (13) 

with N r :=J2 ier f}U 

The correct expression for the partial trace follows from its 
defining property that Tr(^4_B) = Tr(Tr r (^4)B) for all oper- 
ators B that have no support on modes r. Hence, let us con- 
sider such an operator B : J nUr — > J- m ur with no support on 
r, i.e., fiB± = ±B±fi Vi Sr . Let m, n, r be orderings for the 
modes in m, n, and r. Further let a — m © r and b = n © r. 
The operator's matrix elements obey 

a (mr'\B\nr) b 

= a (0\K'F™B(P:)HF^\0) b 

= (-i) fWf " m (0| r <0|Fr^5(F;)t(F n ")t|0) r | ) n 

= <W(-l) f( ™ +f ° m(m\B+ + (-l) F 5_|n) n 

= 8 rr i m {m\B\n) n . (14) 

Requiring that 

Tr(iB) = b (nr\A\mr) a m (m\B\n) n = Tr((Tr r A)B), 
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Figure 4: (Color online) Listing of all (contraction) operations that are needed to evaluate a FOC: (a) partial multiplication B •„ A, (b) partial 
trace Try A, (c) partial contraction B o A = Try £> •„ A, and (d) partial projection. The latter two operations are not primitives, but are rather 
applications of partial multiplication and trace. For numerical purposes, it is however useful to implement them. In all cases, the lower operator 
is defined to come first in the operator ordering. 



is true for all operators B with the properties stated above, 
leads to the conclusion that the partial trace for the modes r is 
simply given by the expression 

Tr r A = ^2\n) nb (nr\A\mr) am (m\, (15) 

r 

Hence, assuming we have the operator in the representation 
Jb, a (A), the resulting operator Try A : J- m — > J- n is in the 
occupation number representation 

(n\J n , m (Tr r A)\m) = (nr\J b ^(A)\mr) . (16) 

Please note that we have chosen the orderings of the modes 
such that, in Eq. ( |14) , two sign factors compensate - that of a 
mode reordering with one from commuting FJ and the oper- 
ator B. A sign factor (— iy( m + n ) would occur in the expres- 
sions for the partial trace, had we swapped the order of m (n) 
and r in the ordering of the incoming (outgoing) modes, i.e., 
a = r©m(b = t0n) instead of our choice here. For such a 
case, the preparative mode reordering would take account of 
the sign factor and then, having realized the preconditions of 
it, one would apply rule ( [To} . 

D. Contraction of some outgoing modes of A with the 
corresponding incoming modes of B and vice versa 

Combining partial multiplication (JTOj with partial trace 
( fT5| l we obtain a general partial contraction, namely, that of 
some outgoing modes n of operator A with the correspond- 
ing incoming modes of B and, simultaneously, contraction of 
some outgoing modes r of B with the corresponding incom- 
ing modes of A. This corresponds to the partial contraction 
depicted in Fig.|4j;. 

Let A : F mUr -> F nUp and B : F nUq -> T kUr , i-e., the 
operators outgoing/incoming supports overlap in the modes n 
and r. Let m, n, r, p, q, t be orderings for the modes in rn, 
n, r, p, q, and fc. Assuming we have the two operators in 



representations A = J n e P ,mer(-4) and B = J £ffirineq (B), 
with a = t ffip and b = m © q, the resulting operator C : 

J~ mUq J~kUp IS 

C = Tr r B- n A 

_ (_i)M+f(p+?) . \ kp)a (kr\B\nq)(np\A\mr) b (mq\, 

i.e., 

{kp\J a>b (C)\mq) 

= (_i)w+f(p+?) . (kr\B\nq){np\A\mr). (17) 

In the following, B o A denotes a (partial) contraction of all 
common incoming/outgoing modes of A with corresponding 
outgoing/incoming modes of B according to Eq. ( [17} . 

E. Partial projection 

The partial projection for an operator is depicted in Fig.[4}l. 
Let A : T m — > J>un- Let r, m, n be orderings for the modes 
in r, m, and n. Further let a = r0n. After projection of modes 
r onto a basis state ({fii} ier eigenstate) \r') t — (F£ y\0) r , 
the resulting operator A' : F m — > T n is 

k = ,(r'\-\rn) a (rn\J atm (A)\m) m (m\ 

— \ n )n (m|, (18) 

i.e., 

(n\J n , m (A')\m) = {r'n\J a , m (A)\m). (19) 

A sign factor (— l) r ™ would occur, if we would swap the or- 
der of modes r and n in the order a of the outgoing modes. 

IV. OPERATOR ORDER AND CONTRACTION 
SEQUENCE 

In Sec. |IIB| the value of the FOC was defined as the 
value resulting from executing the contractions of the con- 
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C 



n 

k t 

J i T 



B 



T 



Figure 5: The most general FOC with three operators. To verify that 
the contraction of operators as given by rule {17} is associative, one 
needs to compare the results of C o (B o A) and (CoB)o A. Both 
do indeed agree. 



stituting operators Ai with respect to a certain operator or- 
der, An o . . . o A-x o Ay. This definition is only sufficient 
if the contraction ([17} of operators, as depicted in Fig. |4j;, is 
indeed associative. For the most general FOC of three oper- 
ators A : J- aU bUc —> FdUeUf' B : J-fugUh — > -^cUfcUn. an d 

one finds indeed (see Fig. [5} 



C o (B o A) = (C o B) o A, 



(20) 



confirming the consistency of the contraction rule fP7} . 

Numerically it may be more efficient to execute for exam- 
ple first the contraction between A% and A3 and contract the 
result with A2 afterwards. To be able to choose an arbitrary 
sequence for the contractions as is possible for the correspond- 
ing QUOCs, we need to be able to change the ordering of the 
operators without changing the value of the FOC. In the el- 
ementary contractions, the ordering of the affected operators 
matters, i.e., for two operators A : F m ur -^nUp an d B : 
FnUq ~~ * ^itur. we have in general Tr r B •„ A ^ Tr„ A > B. 
However, if each of the two operators is either parity preserv- 
ing (s = 0) or parity changing (s = 1), we find the simple 
relation 



at the corresponding contraction arc as exemplified in Fig. [6] 
The numerical overhead for keeping track of those signs is 
marginal. 

In the following, operators A : T m — > T n that are either 
fermion number parity preserving or changing, 



(-1) N «A = ±A(-1) 



N„ 



(22) 



are called parity-symmetric. Also FOCs that contain only 
parity-symmetric operators are called parity-symmetric. 

Using the above result, it is possible to do the operator con- 
tractions of a parity-symmetric FOC in an arbitrary sequence. 
One starts with the predefined operator order. To execute the 
contraction of two (arbitrary) operators of the FOC: 

• Apply rule ( f2"Tj ) to bring the two operators into direct 
neighborhood in the operator order, keeping track of the 
resulting sign factors for the contraction arcs and of the 
global sign, 

• apply mode swapping operators as described in 
Sec. |III A| to bring the occupation number representa- 
tions of the two operators into accord with the precon- 
dition of the general contraction rule ( [17} , and 

• replace the two operators by their contraction according 
to the rule {FT}. 

Consequently, the contraction of a FOC can be done effi- 
ciently - with the same sequence of partial contractions as for 
a corresponding qudit operator circuit. No non-local Jordan- 
Wigner transformations occur. Marginal computational over- 
heads result from keeping track of certain sign factors when 
doing contractions in a sequence that deviates from the or- 
dering of the circuit's operators and reordering of modes for 
incoming or outgoing modes of single operators, directly be- 
fore a partial multiplication, trace etc. that they are affected 
by. 

The operator order is part of the definition of a FOC. For 
the example of the fermionic MERA it can be chosen to agree 
with the physical interpretation as consecutive renormaliza- 
tion steps; i.e., the operator order is increasing with the renor- 
malization number. As all unitaries (isometries) of a particular 
renormalization stage commute, the ordering among those can 
be chosen arbitrarily. In Sec. VII a useful operator ordering 



Tr r B •„ A — ( 1) A B Tr„[(P n •„ A > P r ) > B], (21) f or fermionic PEPS is presented 



where P n : T n -> F n with n (n'\P n \n) n 



In 



W(-i) 

the more compact notation this reads B o A = (—l) SASB P n o 
A o P r o B. In an implementation, instead of inserting the P n 
in this fashion as operators or applying them directly to A or 
B, more efficiently, one may introduce a binary counter (with 
initial state 0) for each contraction arc - in this case, for the 
contraction with respect to modes n. Whenever a factor P n 
arises when swapping the order of operators that have both 
support on n, the state of the binary counter is inverted. Once, 
the contraction with respect to modes n is executed, one in- 
serts the factor (—1)" in the corresponding expression, if the 
state of the counter is 1 . In the graphical representation of a 
FOC, we denote the state 1 of the counter by a minus sign 



V. COMPUTATIONAL COSTS AND LOCALITY 

Given a contraction sequence for a qudit operator circuit 
(QUOC), the same sequence can be used for a correspond- 
ing parity-symmetric FOC (for which all qudit operators are 
replaced by parity-symmetric fermionic operators of identi- 
cal dimension). There is hence no memory or computational 
overhead per se. For the elementary contraction operations 
stated in Sec. Ill a certain ordering of the modes was being 
assumed, prior to the operation. If one uses the contraction 
operations as stated there, one gets a marginal overhead from 
the corresponding preparative mode reorderings; Sec. Ill A 
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Figure 6: (Color online) To allow for arbitrary contraction sequences, one needs to be able to change the operator ordering. In the diagrams, the 
operator order is defined to increase from the bottom to the top. If each operator is parity-symmetric (either preserves or changes the fermion 
number parity; s = or s = 1), swapping of operators can be done and the resulting sign factors taken account of efficiently, (a) The generic 
rule l |2 1 [ i for swapping two operators that are neighbors in the ordering, (b) Identities for the most generic FOC with three operators, the same 
as in Fig. [5] depicted in a slightly different fashion. A minus sign at the contraction arc for a mode set n indicates that a sign factor (—1)" is 
to be inserted in the contraction formula (see text). 



The number of numerical operations needed for a reorder- 
ing is proportional to the size of the operator matrix: every 
reordering can be achieved by a sequence of swaps of con- 
secutive modes. The product of appropriate swaps yields a 
reordering operator that is sparse with exactly one entry ±1 
in each row and column. To apply such an operator to either 
side of A : T m — > T n , requires only XmXn operations, where 
Xm and Xn are the dimensions of the (possibly reduced) in- 
coming and outgoing Hilbert spaces. Every contraction of the 
operator, except for partial traces or projections, would how- 
ever already require a larger number of numerical operations. 
The computational overhead is hence marginal. There is no 
overhead in memory requirements. 

Further, all considerations about locality, hence, carry over 
directly from those of the known QUOCs (for instance the 
qudit MERA) to the corresponding FOC (e.g., the fermionic 
MERA). In the calculation of local expectation values w.r.t. a 
MERA, only operators inside a causal cone of the observable 
enter the actual calculation (all others cancel). That Jordan- 
Wigner strings outside the causal cone can be avoided for the 
fermionic MERA has already been shown by an alternative 
approach in Ref. |3Q| . see also Ref. |29|. 



B. Reversing contraction arcs 

For algorithms operating on FOCs, as for example the one 

it is sometimes use- 



for fermionic PEPS presented in Sec. VII 



ful to reverse contraction arcs, i.e., to change outgoing modes 
of one operator to incoming modes and vice versa at the op- 
erators it is contracted with; see Fig. [7^. Let A : J-'mur ~^ 
J 7 , 



nUsUp 



and B : J- nUsUq —> J-kur, i-e., the operators outgo- 
ing/incoming supports overlap in the modes n, r and s. Let 
m, n, r, s, p, q, t be orderings for the modes in m, n, r, s, 
p, q, and k. For reversing the arc corresponding to modes n, 
i.e., changing the modes n to be incoming (outgoing) at oper- 
ator A (B), the relations between A and B and the resulting 
operators (as depicted in Fig. [7^) are 



tffir (fcr|B|nsq) nffi5ffiq = (-1)'' 
nme ® p {nsp\A\mr) m@t = (-1)' 
such that B o A = B' o A'. 



n®x(knr\B'\sq) s(Sq , 
5 p (sp\A'\mnr 



mffitiffiv 1 



C. Singular value decomposition and truncation 



VI. FURTHER OPERATIONS ON FOCS 
A. Hermitian conjugation 

The Hermitian conjugate of a FOC is simply given by 

(A N o . . . o li)t = A{ o . . . o A\. (23) 



The operator order is reversed and one has to take the Her- 
mitian conjugate of each fermionic operator in the circuit. In 
the representation as a directed graph, all arcs are reversed. 
The Hermitian conjugate is for example of interest when cal- 
culating expectation values with respect to a (pure) FOC state. 
Fig. [8^ shows it for the example of a fermionic PEPS. 



It is possible to decompose an operator A : J- m un —> ^uuv 
by singular value decomposition with respect to arbitrary 
splittings of the incoming and outgoing modes. The resulting 
circuits can be chosen to be of the form C0B01C0A0B, 
where A : T z T x = \z\) is a diagonal operator en- 
coding the singular values; see Figs.[7]3-[7}l. This also allows 
for truncation of modes (or the reduction of Hilbert space di- 
mensions): Contract two operators C : F m ux J~ u and 
B : T n — > J- X uv, as in Fig. [7c to obtain an operator A, apply 
the singular value decomposition to it and truncate (some of 
the smallest) singular values, to obtain an approximation of 
C o B where the dimension of the retained Hilbert space for 
the modes in x has been reduced. 

Let m, n, u, 0, y, and 3 be orderings of modes in m, n, u, 
v, x, and z. The contraction of the FOC C o B, as depicted in 
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(b) 



A 



(c) 



C 



B 



m n 



(d) u v 
i_ 



(7 



z 
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Figure 7: (Color online) In all subplots, the operator order is denned 
to increase from the bottom to the top. (a) It is possible to reverse 
contraction arcs. The resulting operators can be expressed in terms 
of matrix elements of the original operators; see Sec. |VIB| Reversing 
the arc for modes n yields the sign factor (— (b-d) It is pos- 
sible to decompose operators (A) by singular value decomposition, 
resulting in circuits of the form (c) or (d). This also allows for the 
reduction of retained Hilbert space dimensions: Contract operators 
B and C to obtain an operator A, apply the singular value decom- 
position to it and truncate (some of the smallest) singular values, to 
obtain an approximation of C o B. Reversing contraction arcs and 
truncation of Hilbert spaces via singular value decomposition are for 
example employed in the contraction algorithm for fermionic PEPS 
in Sec. EH] 



Fig. [7J: yields 

CoB = (-l)™\uv) um u (u\C\xn) f(Sn 

x (xv \B\n) n n (g m (nm\. (24) 

With the occupation number representation A := 
•Aieo.nemC^) of A, we can hence decompose the oper- 
ator by applying the singular value decomposition to the 
matrix A defined by 

(uv\A\nm) := (-1)™ (uv\A\nm), (25) 
A = UAV, (26) 

where U and V are unitary and A is the diagonal matrix of 
singular values. The operators of the resulting circuit CoB 
can then be chosen as (0 < a < 1) 



Ju, f ®m{C) = UA a , J mo , n (B) = ^~ a V. 



(27) 



When the singular values are to be separated into a third op- 
erator A as depicted in Fig. [7J1, the operators of the resulting 



circuit C o A o B are given by 

Ju***{p) = U, J M (A)=A, J 3 ®o,n(B) = V. (28) 

Reduction of Hilbert space dimensions {truncation) via sin- 
gular value decomposition is for example employed in the 
algorithm for evaluating expectation values with respect to 



fermionic PEPS in an approximative fashion; see Sec. VII 



VII. FERMIONIC PEPS 

The FOC framework incorporates a fermionic version of 
the class of qudit states called tensor product ansatze ifTBT - 
[T8l or projected entangled pair states (PEPS) lfl9l . In Ref. 
||3T|| it was suggested to obtain fermionic PEPS by applying 
fermionic parity-symmetric (projection) operators to a ten- 
sor product of maximally entangled pair states. The detour 
over maximally entangled states is not necessary (but also 
not harmful); as depicted on the left hand side of Fig. [8^, 
a fermionic PEPS on a square lattice, equivalently, can be 
defined by assigning to each lattice site (x, y) (away from 
the boundaries) a parity-symmetric fermionic operator A : 
J oUm — > Tbunus where a and m are sets of incoming modes 
from operators on neighboring sites (x + 1, y) and (x, y + 1), 
and b and n are outgoing modes to operators on sites (x—l,y) 
and (x,y—l). The set of modes s composes the local physical 
Hilbert space of site (x,y). In the FOC framework, the gener- 
alization to more complicated or higher-dimensional lattices 
is straightforward. The choice of the direction of the con- 
traction arcs is (an arbitrary) part of the definition of the state 
and can also be changed later as described in Sec. 



VI B To 



complete the definition of the fermionic PEPS one needs to 
specify an (initial) operator order. An example is given on 
the left hand side of Fig. [8^, where the gray line below the 
lattice indicates the lexicographic order with respect to lattice 
coordinates (— x, y). 

In Ref. ||3T1 it was described how the FOC of a fermionic 
PEPS can be mapped to a QUOC by choosing a fixed ordering 
of all modes. This was achieved with one additional bond per 
horizontal contraction arc (i.e., a factor of four in the number 
of degrees of freedom per site) and a correspondingly reduced 
computational efficiency (a factor of several powers of four) 
for the evaluation of expectation values, calculation of ground 
states etc. 

The approach presented here is an alternative one, empha- 
sizing that the mapping to a QUOC (with a fixed mode or- 
der) is not necessary. All manipulations and contractions 
on fermionic PEPS can be done according to the rules de- 



scribed in Sees. Ill IV and VI In that case, compared to the 
same operations on a corresponding qudit PEPS (replacing the 
fermionic operators with qudit operators of identical dimen- 
sions), only marginal computational overheads arise. 

Fig.|8]shows graphically how the FOC for the evaluation of 
a local expectation value (tp\0\ip) can be constructed. For the 
bra vector (dual vector) (ip\, operator order and contraction 
lines reverse as a side effect of taking the Hermitian conju- 



gate; Sec. VI A For later convenience this is reverted by ap- 
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Figure 8: (Color online) (a) A fermionic PEPS can be constructed as a FOC, where fermionic operators are assigned to each lattice site. As 
chosen here for a square lattice, each operator has two sets of incoming modes from operators on neighboring sites and two outgoing sets 
of modes to operators of the remaining nearest neighbors. One outgoing set of modes corresponds to the physical site Hilbert space. The 
Hermitian conjugate of the circuit An ° • ■ • ° Ai is A\ o • ■ • o A N . All contraction arcs and the operator order (gray line below/above the 
circuit) are reversed. This side effect can be reverted (without changing the value of the FOC) by applying Eq. l |21| > and the rule derived in 
Sec. |VIB] with only a marginal computational overhead, (b) To evaluate a local expectation value, the FOCs for bra, local observable, and 
ket have to be composed. The operator order can again be changed for later convenience - in this case no additional sign factors occur, as all 
swapped operators have no common contraction arcs. For the definition of the objects on the right hand side, see also Fig. [9^. 
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Figure 9: (Color online) (a) Definition of the objects on the right hand side of Fig.JSJj - here, in particular, for the site where the local observable 
acts nontrivially. (b) The FOC for the evaluation of a local observable is contracted by considering the first row of the FOC as a fermionic 
state \xi) and applying the other rows as operators to it \\ y ) = T y \x y -i)- Doing this in an exact manner, the number of degrees of freedom 
per site for the states \x v ) would in general increase exponentially with y. One can decrease them during the algorithm for the case of a finite 
(infinite, translationally invariant) lattice by applying the DMRG (iTEBD) algorithm. 
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plying Eq. pTj ) and the rule derived in Sec. VI B After com- 
posing bra, observable, and ket, the operator order can again 
be changed conveniently, this time without any sign factors 
occurring, as all swapped operators share no common contrac- 
tion arcs; Fig. [8(3. As in the qudit case 1191 . the contraction of 
the resulting circuit can be executed row by row, i.e., by treat- 
ing the lowest row as a one-dimensional fermionic state \xi) 
to which the operators of the following row T y (row trans- 
fer matrix) are applied; \x y ) — T y \x y -i)- No additional 
sign factors occur due to operator reorderings (Fig. 15), but 



only due to mode reorderings (Sec. Ill Ai before contractions 



(marginal overhead). An essential aspect of PEPS algorithms 
is that contractions, e.g., for the evaluation of expectation val- 
ues, cannot be executed exactly, as the stepwise application 
of the row transfer matrices would in general lead to an expo- 
nential growth in the number of modes per site for \Xy}- As 
suggested in Ref. [ 19], this can be circumvented by applying a 
variant of the density-matrix renormalization-group (DMRG) 
algorithm fl4j[5) to each state \Xy)> before executing the con- 
tractions to the next row. The only purpose of the DMRG 
procedure is here to reduce the number of degrees of freedom 
in each step to a manageable number, and hence, do contrac- 
tions in an approximative fashion. The essential operation is 
to do Schmidt decompositions of \xy)- This can be done for 
FOCs as described in Sec. |VTCl 

The FOC framework also allows to simulate infinite 
fermionic PEPS. To this purpose, the fermionic PEPS is to 
be defined by repetition of an elementary cell FOC; cmp. to 
Ref. Il33l for the qudit case. The algorithm does not deviate 
substantially from the finite-size case. The biggest difference 
being that, for the reduction of degrees of freedom in states 
\Xy)> one nas to use a translationally invariant formulation of 
the DMRG algorithm, basically the iTEBD algorithm as de- 
scribed in Ref. [38], again based on the ability to do singular 
value decompositions (Sec. |VIC) . With this, one has a trans- 
lation of the algorithms for the calculation of approximative 
ground state or time-evolved qudit (i)PEPS |[T51 l33l to the 
fermionic case without reduction of the computational effi- 
ciency, as those algorithms are based on the ability to contract 
operator circuits just as in our example. 



VIII. DISCUSSION 

In Ref. 11301 it was shown that contractions of fermionic 
unitary circuits with a causal cone (for instance the evalua- 
tion of local observables w.r.t. a MERA) can be done without 
occurrence of any Jordan- Wigner strings outside the causal 
cone. Here, this result was extended in proving that arbitrary 
parity-symmetric fermionic operator circuits can actually be 
contracted with the same computational effort and memory 
requirements as a corresponding QUOC. This remarkable re- 
sult follows from the fact that a given contraction sequence for 
a QUOC can be implemented for a corresponding FOC with 
essentially the same number of computational operations. We 
have presented the required contraction primitives and dis- 
cussed the marginal computational overheads. 

This allows to translate algorithms on QUOCs to corre- 



sponding algorithms on FOCs. For example in the algorithm 
for scale-invariant MERA as studied in Refs. (27l|351|40), the 
super operator simply becomes a fermionic super operator. Its 
iterative application to an observable yields the expectation 
value of the observable in the thermodynamic limit. 

For the special example of the FOC being a MERA, in Ref. 
E51l , first numerical results where presented (postponing a de- 
scription of the algorithm for a later publication). A scheme 
for fermionic PEPS was suggested in Ref. OTI . The suggested 
mapping to a QUOC used there seems numerically less effi- 
cient than the contraction scheme presented here. Instead of 
encoding the fermionic sign factors by increasing tensor di- 
mensions, they can be taken account of during contractions, 
specifically in preparative mode reorderings, and operator or- 
der swaps. The resulting marginal overhead appears smaller. 

It will be interesting to see to what extent variational 
ansatze like fermionic variants of PEPS or MERA, both satis- 
fying entropic area laws BTH431I . will be able to appropriately 
grasp the correlations present in critical fermionic strongly 
correlated models, models that are known to violate such area 
laws logarithmically [44-48]. First numerical results ||25H3TI 
seem promising. It is the hope that the framework discussed in 
this work will help in constructing fermionic variants of vari- 
ational approaches to simulate strongly correlated fermions in 
higher dimensions. 
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Appendix A: Alternative derivation of the rule for partial 
multiplications 



For the same preconditions as in Sec. MB we want to de 



rive the partial multiplication rule - this time by extending the 
supports of operators A and B prior to the contraction. 

Let D : T m — > F n , and let m, n, r be orderings of the 
modes in m, n, and r. Extending the incoming and outgoing 
supports of D by modes r (on which the resulting operator is 
supposed to act trivially), we arrive at D' : J- m ur — > J- nU r 
with 

D = I TJ.r) n t n (n\D\m) m m(Sv( mr \- (Al) 

This is confirmed by Tr r D' = D, according to ( [T5) , and 
f x D' ± =±D' ± f x for dlxtr. 

Using the rule ( |A 1 1 > to extend the supports of the operators 
A : T m — > J-nup an d B : J- n uq —> J~k by modes q and 
p, respectively, defines operators A' : J- m u q —> J'nupuq an d 
B' : J- n upuq — > Fkup- The partial multiplication C = B - n A 
amounts now simply to the usual operator product C = B'-A' . 
Assuming we have the two operators in representations A = 
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Ja,m(A) and B — Jf ib (B) with a = n ffl p and b = n © q (p The result C : J r m ug —> J~kup of the multiplication with 

and q are orderings for the modes in p, and q), orderings Ci = t © p, C2 = m © q is then 

A' = \npq) am a {np\A\m) m mSq {mq\, q _ & ■ A' 

and = (-iflMtep i k \ B \ n ^) ( J l P\ A \ rn ) m<sq{ m( l\ 

= \kp)c x {kp\C\mq) t2 {mq\, 

B' = |fep)tep t(k\B\nq) b bmp (nqp\ 

= (-l) pq \kp) t(Sp (k\B\nq) nfBmq (npq\. coinciding with Eq. (TO), 
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